Leakage Belt Delineation

VM0048-compliant delineation of valid leakage area for the Gola REDD+ Forest Carbon Project project area expansion

Authors

Affiliation

Seamus Murphy

Winrock International

Krysla Lima Fernandes

Winrock International

Rajesh Bista

Winrock International

Meyru Bhanti

Winrock International

Project Area & Jurisdictions

country = sf::read_sf("./data/AOI/liberia_boundary_national.shp") |>sf::st_transform(32629)
counties = sf::st_read("./data/AOI/places_poly_county.shp") |>sf::st_transform(32629)

Reading layer `places_poly_county' from data source 
  `/Users/seamus/repos/rspb-redd-risk-new/data/AOI/places_poly_county.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 16 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -11.50675 ymin: 4.353908 xmax: -7.367323 ymax: 8.551925
Geodetic CRS:  WGS 84

jurisdiction = counties |>dplyr::filter(name=="Grand Cape Mount County"|name=="Gharpolu County")
jurisdiction$name = 'Grand Cape Mount & Gharpolu Counties'

aoi = sf::read_sf("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/AOI/Archive/ProjectArea_CF-Expansion_051525/updated_areas.shp") |>
  sf::st_make_valid() |>
  sf::st_transform("EPSG:32629") |>  
  sf::st_cast("MULTIPOLYGON") |> sf::st_as_sf()# |> dplyr::select("Name")

aoi$area_ha = round(as.numeric(sf::st_area(aoi) * 0.0001, 4))
aoi |> sf::st_drop_geometry() |> janitor::adorn_totals() |> 
  flextable::flextable() |> 
  flextable::fontsize(size=8,part="all") |> 
  flextable::autofit()

Name	Shape_Leng	Shape_Area	area_ha
Lukasu	138,549.68	469,973,056	46,997
Mbarma CF	91,209.87	443,230,136	44,323
Lower Sokpo CF	59,692.55	140,691,204	14,069
Upper Sokpo CF	62,831.33	109,225,715	10,923
Zuie CF	92,221.99	362,360,113	36,236
Gola Forest National Park	182,951.78	890,423,315	89,042
Tonglay CF	76,824.36	195,953,085	19,595
Norman CF	77,870.76	123,489,968	12,349
Foya	351,458.82	1,048,825,566	104,885
Kposo	77,241.08	307,632,539	30,763
Gola Konneh CF	117,719.73	531,018,504	53,102
Gola Konneh CF (ext. S)	56,323.80	108,968,871	10,897
Gola Konneh CF (ext. N)	96,399.13	285,898,112	28,590
Lower Sokpo CF (ext.)	45,714.33	47,712,283	4,771
Maima CF (ext.)	59,344.81	139,163,854	13,916
Zuie CF (ext.)	17,920.40	14,953,753	1,495
Yangaya	46,462.77	58,336,055	5,834
Koninga	71,795.27	317,410,805	31,741
Lobarsu	62,061.63	165,254,549	16,548
Bade	99,147.36	513,666,252	51,367
No-mans land	209,804.70	394,466,117	39,447
Lower Guma	56,222.96	142,192,367	14,219
Central Guma	31,669.60	38,501,602	3,850
Upper Guma	54,605.75	137,757,843	13,771
Hassala	62,847.92	145,701,577	14,570
Buluyeama	47,779.35	100,212,025	10,021
Hembeh	199,395.60	461,835,431	46,184
Total	2,546,067.34	7,694,854,697	769,505

Derive Leakage Area Belt

aoi_union = sf::st_transform(aoi, 32629) |>  sf::st_union() |> sf::st_make_valid()
leakage_buffer = sf::st_buffer(aoi_union, dist = 5500, endCapStyle="ROUND") |>
  sf::st_as_sf()
leakage_buffer = concaveman::concaveman(leakage_buffer, concavity=5) |> 
  sf::st_zm() |> sf::st_make_valid()
leakage_belt_whole = sf::st_buffer(leakage_buffer, dist = 4500, endCapStyle="ROUND") |>
  sf::st_make_valid() 

tmap::tmap_mode("view")
tmap::tm_shape(leakage_belt_whole) + 
  tmap::tm_polygons(col="orange",fill="orange",fill_alpha=0.3,lwd=1)+
  #tmap::tm_add_legend(type="lines",col="orange",labels="Leakage Belt (10km)") +
  tmap::tm_shape(aoi) + tmap::tm_borders(lwd=1, col="white") +
  tmap::tm_basemap("Esri.WorldImagery")

Reading layer `LeakageBelt_10k-Radius_UnFiltered-Whole' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered-Whole.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -11.31102 ymin: 7.024169 xmax: -9.847199 ymax: 8.24916
Geodetic CRS:  WGS 84

leakage_belt    = sf::st_difference(leakage_belt_whole, st_union(st_combine(aoi_union)))
leakage_belt$area_ha = round(as.numeric(sf::st_area(leakage_belt) * 0.0001, 4)) 
leakage_belt = sf::st_intersection(country, leakage_belt)

tmap::tmap_mode("view")
tmap::tm_shape(leakage_belt) + 
  tmap::tm_polygons(col="yellow",fill="yellow",fill_alpha=0.3)+
  #tmap::tm_add_legend(type="lines",col="yellow",labels="Leakage Belt (10km)") +  
  tmap::tm_shape(aoi) + tmap::tm_borders(lwd=0.4, col="white") + 
  tmap::tm_text(text="Name", size=0.5, col="white") +
  tmap::tm_basemap("Esri.WorldImagery")

# save locally
#sf::wt_write(leakage_belt, "OneDrive.../20087 - RSPB Gola Feasibility/Deliverables/
#  Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered.zip") 
#sf::st_write(leakage_belt_whole, "OneDrive.../20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered-Whole.shp")

Reading layer `leakage_belt_10km_liberia_ext' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered/leakage_belt_10km_liberia_ext.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 245117.9 ymin: 776808.3 xmax: 406590.1 ymax: 908414.4
Projected CRS: WGS 84 / UTM zone 29N

Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 245117.9 ymin: 776808.3 xmax: 406590.1 ymax: 908414.4
Projected CRS: WGS 84 / UTM zone 29N
                      Name Shape_Leng Shape_Area                       geometry
1 Leakage Belt 10km Radius   138549.7  469973056 MULTIPOLYGON (((245306.7 79...
  area_ha
1  298838

Derive Leakage Area Masks

Roads Mask

roads_one = sf::st_read("~/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_one.shp")
roads_two = sf::st_read("~/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_two.shp")

# we have simplify mask shapefiles and split them up to shorten computing 
# time & avoid crashing. See option for "harsh" simplification on line 163
roads_one_simplified = roads_one |> sf::st_make_valid() |> sf::st_cast("MULTILINESTRING") |> 
  rmapshaper::ms_simplify(keep=0.5)
roads_two_simplified = roads_two |> sf::st_make_valid() |> sf::st_cast("MULTILINESTRING") |> 
  rmapshaper::ms_simplify(keep=0.5)

# bigger file needs more simplificaiotn
roads_one_simplified_harsh = rmapshaper::ms_simplify(
  roads_one_simplified, keep=0.01) 

# now apply buffer operation, but note this takes time. Its 
# advised processing inputs as much as possible before running
roads_one_buffer = sf::st_buffer(
  roads_one_simplified_harsh, 
  dist = 10000, 
  nQuadSegs = 5,
  endCapStyle="ROUND", 
  joinStyle = "ROUND",
  mitreLimit = 1,
  singleSide = FALSE
  )

roads_two_buffer = sf::st_buffer(
  roads_two_simplified, 
  dist = 10000, 
  nQuadSegs = 5,
  endCapStyle="ROUND", 
  joinStyle = "ROUND",
  mitreLimit = 1,
  singleSide = FALSE
  )

# Combine, dissolve and cast to single feature
roads_mask = sf::st_combine(roads_one_buffer, roads_two_buffer) |>
  sf::st_union() |> sf::st_cast("POLYGON")

# Visual check
tmap::tmap_mode("view")
tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=0) +
  tmap::tm_shape(roads_one_simplified_harsh) + tmap::tm_lines(lwd=2, col="red") +
  tmap::tm_shape(roads_two_simplified) + tmap::tm_lines(lwd=2, col="green") +
  tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=1, col="pink") + 
  #tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_scale_bar(position = c("RIGHT", "BOTTOM"), text.size = .5) + 
  tmap::tm_compass(color.dark = "gray60", text.color = "gray60", position = c("left", "top"))

# Save output to MASKS folder and purge memory
#sf::st_write(roads_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASKS/LeakageMask_Roads_10km-Buffer_051625.shp", delete_dsn=T)

Reading layer `roads_simplified_one' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_one.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 842 features and 16 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 237648.4 ymin: 761962.5 xmax: 400449.6 ymax: 923370.5
Projected CRS: WGS 84 / UTM zone 29N

Reading layer `roads_simplified_two' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_two.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 255 features and 16 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 237648.4 ymin: 762265 xmax: 400334.4 ymax: 923370.5
Projected CRS: WGS 84 / UTM zone 29N

Reading layer `Road_Mask_10km-Proximity_051625' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/Archive/LeakageMask_Roads_10km-Buffer_051625/Road_Mask_10km-Proximity_051625.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 17 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 227660.4 ymin: 752279.5 xmax: 410286 ymax: 933257.1
Projected CRS: WGS 84 / UTM zone 29N

Habitat Mask

# import inputs
wetlands = terra::rast("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Wetlands/GLWD_EPSG32629.tif")
protected_areas = sf::st_read("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Protected Areas/Archive/WDPA_Mar2025_Public_32629_GOLA.shp")

Reading layer `WDPA_Mar2025_Public_32629_GOLA' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Protected Areas/Archive/WDPA_Mar2025_Public_32629_GOLA.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 6 features and 30 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 240189.7 ymin: 751293.6 xmax: 411102.7 ymax: 895314.1
Projected CRS: WGS 84 / UTM zone 29N

leakage_belt_crop = sf::st_transform(sf::st_as_sf(leakage_belt_whole), 32629) |> terra::vect()
wetlands          = terra::crop(wetlands, leakage_belt_crop, mask=T)

# tidy labeling
code_dict_2 <- data.frame(
  id = c(1, 4, 7, 10, 12, 14, 15, 18, 20, 21, 26, 31),
  label = c(
    "Freshwater lake",                              # 1
    "Large river",                                  # 4
    "Small streams",                                # 7
    "Riverine, regularly flooded, forested",        # 10
    "Riverine, seasonally flooded, forested",       # 12
    "Riverine, seasonally saturated, forested",     # 14
    "Riverine, seasonally saturated, non-forested", # 15
    "Palustrine, seasonally saturated, forested",   # 18
    "Ephemeral, forested",                          # 20
    "Ephemeral, non-forested",                      # 21
    "Tropical peatland, forested",                  # 26
    "Other coastal wetland"                         # 31
  ))

levels(wetlands) <- code_dict_2
wetlands[wetlands == 0] <- NA

# derive wetland mask
wetlands_mask <- wetlands
wetland_classes <- c(1, 4, 7, 10, 12, 14, 15, 18, 20, 21, 26, 31)
terra::values(wetlands_mask) <- ifelse(terra::values(wetlands) %in% wetland_classes, 1, NA)

# save locally for faster computing
#raster::writeRaster(wetlands_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_Wetland-GLWD_051625.tif", overwrite=T)
#sf::st_write(protected_areas, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_ProtectedAreas_WDPA_051625.shp", delete_dsn=T)

Slope Mask

# skipping these operations here (est. time 12 mins)
DEM = terra::rast("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/DEM/DEM_EPSG32629.tif") 

# derive slope percentage from degree 
slope_degrees = terra::terrain(DEM, v="slope", unit="degrees")
slope_percent = tan(slope_degrees * (pi / 180)) * 100
slope_percent = terra::clamp(slope_percent, 0, 100) 
slope_invalid = slope_percent > 10
slope_invalid[slope_invalid == 0] <- NA
slope_mask = terra::as.polygons(slope_invalid, dissolve=T)|>sf::st_as_sf()|>sf::st_union()

# save locally & reload
sf::st_write(slope_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_Slope10%-Invalid_051625.zip", delete_dsn=T)
slope_mask = sf::st_read("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/Archive/LeakageMask_Slope10%-Invalid_051625/slope_poly_simplified.shp")

Visual Check

# Visual check
tmap::tmap_mode("view")
tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=0) +
  tmap::tm_shape(wetlands) + tmap::tm_raster(col.legend = tm_legend("Wetlands (GLWD")) +
  tmap::tm_shape(leakage_belt) + tmap::tm_polygons(col="yellow",fill="yellow",fill_alpha=0.4, lwd=1.5)+ 
  tmap::tm_shape(protected_areas) + tmap::tm_polygons(fill="ORIG_NAME", fill.legend = tm_legend("Protected Areas (WDPA)")) +
  tmap::tm_shape(aoi) + tmap::tm_borders(lwd=1.5, col="red") + 
  tmap::tm_text(text="Name", size=0.3, col="black") +
  tmap::tm_shape(roads_mask) + tmap::tm_borders(col="purple") +
  #tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = .5) + 
  tmap::tm_compass(color.dark = "gray60", text.color = "gray60", position = c("left", "top"))

Apply Leakage Area Masks

# clip
roads_leakage_mask = sf::st_intersection(leakage_belt, roads_mask)
slope_leakage_mask = sf::st_intersection(leakage_belt, slope_mask)
wetlands_leakage_mask = sf::st_intersection(leakage_belt, wetlands_mask)
protected_areas_leakage_mask = sf::st_intersection(leakage_belt, protected_areas)

# merge. Btw actually makes better sense to keep these seperate. 
# They are easier to operate seperately due to their size and linework.
leakage_mask_a = sf::st_union(roads_leakage, slope_leakage_mask)
leakage_mask_b = sf::st_union(leakage_mask_a, wetlands_leakage_mask)
leakage_mask_c = sf::st_union(leakage_mask_b, protected_areas_leakage_mask)

# save
sf::st_write(leakage_mask_c, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Leakage Masks/", delete_dsn=T)

Tally Features

road_count_whole = sf::st_intersection(road, leakage_belt_whole)
road_count = sf::st_intersection(road, leakage_belt)
road_length_whole = sum(sf::st_length(road_count_whole)) + sum(sf::st_length(road_count_whole))
road_length = sum(sf::st_length(road_count)) + sum(sf::st_length(road_count))
road_length_whole
road_length

waterways_count_whole = sf::st_intersection(waterways, leakage_belt_whole)
waterways_count = sf::st_intersection(waterways, leakage_belt)
waterways_length_whole = sum(sf::st_length(waterways_count_whole))
waterways_length = sum(sf::st_length(waterways_count))
waterways_length_whole
waterways_length

places_count_whole = sf::st_intersection(places, leakage_belt_whole)
places_count = sf::st_intersection(places, leakage_belt)
places_count_whole
places_count